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Abstract.  The  Lagrangian  prediction  skill  (model  ability 
to  reproduce  Lagrangian  drifter  trajectories)  of  the  now¬ 
cast/forecast  system  developed  for  the  Gulf  of  Mexico  at  the 
University  of  Colorado  at  Boulder  is  examined  through  com¬ 
parison  with  real  drifter  observations.  Model  prediction  error 
(MPE),  singular  values  (SVs)  and  irreversible-skill  time  (IT) 
are  used  as  quantitative  measures  of  the  examination.  Di¬ 
vergent  (poloidal)  and  nondivergent  (toroidal)  components  of 
the  circulation  attractor  at  50  m  depth  are  analyzed  and  com¬ 
pared  with  the  Lagrangian  drifter  buoy  data  using  the  empir¬ 
ical  orthogonal  function  (EOF)  decomposition  and  the  mea¬ 
sures,  respectively.  Irregular  (probably,  chaotic)  dynamics  of 
the  circulation  attractor  reproduced  by  the  nowcast/forecast 
system  is  analyzed  through  Lyapunov  dimension,  global  en¬ 
tropies,  toroidal  and  poloidal  kinetic  energies.  The  results  al¬ 
low  assuming  exponential  growth  of  prediction  error  on  the 
attractor.  On  the  other  hand,  the  g-th  moment  of  MPE  grows 
by  the  power  law  with  exponent  of  3 q /A.  The  probability 
density  function  (PDF)  of  MPE  has  a  symmetrical  but  non- 
Gaussian  shape  for  both  the  short  and  long  prediction  times 
and  for  spatial  scales  ranging  from  20  km  to  300  km.  The 
phenomenological  model  of  MPE  based  on  a  diffusion-like 
equation  is  developed.  The  PDF  of  IT  is  non-symmetric  with 
a  long  tail  stretched  towards  large  ITs.  The  power  decay  of 
the  tail  was  faster  than  2  for  long  prediction  times. 


1  Introduction 

During  the  World  Ocean  Circulation  Experiment  (WOCE), 
the  ocean  velocity  observation  has  been  significantly  ad¬ 
vanced  with  extensive  spatial  and  temporal  coverage  us¬ 
ing  near-surface  Lagrangian  drifters,  RAFOS  floats,  and 
Autonomous  Lagrangian  Circulation  Explorers  (ALACEs). 
Trajectories  of  these  quasi-Lagrangian  drifters  reflect  the 
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whole  spectrum  of  ocean  motions,  including  meso-  and 
submesoscale  eddies,  various  waves,  inertial  and  semi¬ 
diurnal  motions,  and  provide  invaluable  resources  to  estimate 
the  nowcast/forecast  skill  of  high-resolution  regional  ocean 
models  in  terms  of  model  ability  to  simulate  from  synoptic 
to  submesoscale  eddy  movements. 

However,  direct  comparison  between  model  and  quasi- 
Lagrangian/Lagrangian  observations  is  difficult  since  nei¬ 
ther  the  dynamics  of  numerical  models  nor  the  model  input 
(bathymetry,  external  forcing  and  subscale  parameteriztions) 
are  identical  to  the  reality.  It  needs  to  be  determined  if  the 
model-data  difference  comes  from  deficiency  in  modeling 
ocean  physics,  from  some  unessential  imperfection,  or  from 
unrepresentative  data  (Davis  et  al.,  1996).  It  is  clear  that  the 
high-resolution  ocean  model  and  the  Lagrangian  drifter  data 
should  be  compared  only  in  the  statistical  sense. 

The  commonly  used  methods  for  model-data  comparison 
are  listed  as  follows.  The  first  approach  estimates  the  mean 
pseudo-Eulerian  circulation  pattern  and/or  Eulerian  statistics 
from  Lagrangian  trajectories  (e.g.  Figueroa  and  Olson,  1994; 
Davis  et  al.,  1996;  Acero-Schertzer  et  al.,  1997;  Stutzer  and 
Rrauss,  1998;  Garraffo  et  al.,  2001b).  The  second  approach 
is  to  compute  the  Lagrangian  statistics  for  the  real  drifters 
and  the  modeled  synthetic  particles  and  then  to  compare 
these  statistics  by  way  of  statistical  tests  (e.g.  Garraffo  et  al., 
2001a;  McClean  et  al.,  2002).  The  third  approach  focuses 
on  comparison  of  attractors  reproduced  by  a  model  and  de¬ 
tected  from  the  Lagrangian  data  (Chu  et  al.,  2002a).  The 
EOF  technique  is  the  mathematical  tool  of  this  approach. 
The  fourth  approach  which  will  be  applied  here,  is  to  esti¬ 
mate  Lagrangian  prediciton  skill  through  statistics  of  non- 
asymptotic  indicators  (measures),  such  as  the  model  predic¬ 
tion  error,  the  finite  scale  Lyapunov  exponent  (Lacorata  et 
al.,  2001),  the  singular  vectors  (Lorenz,  1965),  stable  and 
unstable  manifolds  (Wiggins,  1992;  Kuznetsov  et  al.,  2002), 
the  irreversible-skill  time  (Ivanov  et  al.,  1994;  Chu  et  al., 
2002c)  and  others.  Early  the  terms  “time  of  predictabil¬ 
ity”  and  “valid  prediction  period”  were  used  instead  of  “the 
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irreversible-skill  time”,  hereinafter  IT,  (Chu  et  al.,  2002a,b). 
In  our  opinion  the  new  name  “IT”,  better  corresponds  to  the 
physical  nature  of  a  model  prediction  skill. 

The  primary  goals  of  the  present  paper  are: 

(i)  to  understand  what  quantitative  measures  may  be  ef¬ 
fectively  applied  for  the  analysis  of  model  ability  to 
predict  the  real  drifter  dynamics.  Such  an  ability  was 
called  “Lagrangian  model  predictability”  by  Mariano  et 
al.  (2002).  Three  local  (non-asymptotic)  measures:  the 
MPE,  SVs  and  IT  will  be  applied, 

(ii)  to  study  Lagrangian  predictability  of  a  high-resolution 
regional  model  through  non-asymptotic  criteria.  We 
will  examine  the  Gulf  of  Mexico  model  with  1/12° 
horizontal  resolution  (Kantha  et  al.,  1999;  Kantha  and 
Clayson,  2000), 

(iii)  to  find  the  fundamental  statistics  of  model  pediciton 
skill  and  parameterize  an  evolution  law  for  the  mean 
prediction  error.  That  in  general  allow  to  detect  when 
model  Lagrangian  predictability  breaks  down. 

All  these  problems  are  in  focus  of  the  modern  predictabil¬ 
ity  studies  of  the  practical  oceanography  and  meteorology 
(Robinson  et  al.,  1999;  Palmer,  2000).  The  knowledge  of  La¬ 
grangian  predictability  of  high-resolution  models  is  very  im¬ 
portant  for  the  numerous  oceanographic  and  ecological  ap¬ 
plications  because  tracer  dispersion  and  pollutant  spread  are 
most  easily  expressed  in  terms  of  Lagrangian  velocity  corre¬ 
lations. 

The  knowledge  of  a  statistical  law  governing  the  average 
growth  of  prediction  errors,  which  initially  have  a  given  or 
even  the  zero  value,  would  allow  us  to  predict  the  expected 
model  skill  and  extending  over  a  given  period  of  time  in  the 
future. 

Note  that  the  high-resolution  oceanographic  modeling  is 
supported  by  the  argument  that  the  models  of  high-resolution 
reproduce  the  mesoscale  eddy  dynamics  more  correctly  ver¬ 
sus  those  of  eddy-permitting  resolution  (Smith  et  al.,  1992). 
Statistics  of  3-D  mesoscale  eddies  dynamics  reproduced 
by  the  high-resolution  models  with  horizontal  resolution  of 
5-10  km  can  be  assessed  only  through  Lagrangian  drifter  ob¬ 
servations. 

The  outline  of  this  paper  is  listed  as  follows.  Section  2 
discusses  the  measures  of  Lagrangian  predictability,  such  as 
the  MPE,  SVs  and  IT.  An  original  approach  developed  in 
the  present  paper  for  computations  of  the  singular  values  and 
singular  vectors  is  adressed  in  Appendix  A. 

The  spatio-temporal  structure  of  the  model  attractor  is  ana¬ 
lyzed  through  Moffat-Zeldovich  decomposition  (MZD),  em¬ 
pirical  orthogonal  functions  and  global  entropy.  This  is  dis¬ 
cussed  in  Sect.  3. 

Using  MZD  two  subspaces  generated  by  the  divergent 
(poloidal)  and  nondivergent  (toroidal)  EOFs,  respectively, 
are  introduced.  Knowledge  of  EOF  decomposition  for  the 
circulation  allows  estimating  Lyapunov  dimension  of  the  cir¬ 
culation  attractor.  Section  4  contains  the  results  of  such  an 
analysis. 


Statistics  of  the  model  prediction  error  and  the  IT  are  dis¬ 
cussed  in  Sect.  5.  We  demonstrate  that  the  process  of  dis¬ 
placement  with  time  of  the  real  drifters  and  modeled  particles 
is  described  by  non-Gaussian  statistics  for  both  small  and 
large  times.  The  growth  of  MPE  holds  the  power  law.  This  is 
a  signature  of  long-term  correlations  between  the  prediction 
and  reality.  A  simple  empirical  model  is  proposed  for  de¬ 
scription  of  MPE  growth.  We  pointed  out  that  linear-tangent 
models  of  MPE  evolution  are  not  sufficient  for  the  analysis 
of  MPE  in  the  Gulf  of  Mexico  circulation  model. 

Section  6  develops  and  illustrates  a  new  approach  for  the 
reconstruction  of  circulation  from  rare  Lagrangian  observa¬ 
tions.  The  last  section  is  a  summary  and  discussion. 


2  Measures  of  Lagrangian  predictability 

2.1  Local  Lagrangian  predictability 

We  may  speculate  three  possible  scenarios  of  Lagrangian 
predictability.  First,  the  model  reproduces  the  pattern  of 
the  real  circulation  attractor  including  its  small-scale  details 
and  predicts  a  drifter  behavior  for  long  times  (the  global  La¬ 
grangian  predictability).  Second,  in  general  the  circulation 
attractor  is  correctly  reproduced  but  not  all  the  small-scale 
details  of  the  circulation  topology  are  resolved.  Here,  the 
predictability  of  the  drifter  dynamics  exists  only  for  the  short 
and  intermediate  time  intervals  (the  local  Lagrangian  pre¬ 
dictability).  Third,  the  predicted  attractor  pattern  differs  from 
the  observed  one  even  in  large  scale  details  and  in  principle 
the  model  can  approximate  drifter  trajectories  only  within 
very  short  time  intervals  (the  partial  Lagrangian  predictabil¬ 
ity).  Identification  of  the  scenario  of  model  Lagrangian  pre¬ 
dictability  is  an  important  task  in  model  verification. 

In  the  present  study  we  a  priori  assume  that  the  modeled 
and  observed  attractors  are  similar,  i.e.  at  least  their  Lya¬ 
punov  dimensions  coincide.  However,  we  can  doubt  that 
a  model  correctly  reproduces  all  attractor  details  including 
vertical  motions.  Thus,  the  second  scenario  of  Lagrangian 
predictability  is  expected  to  appear  in  our  computations. 

Briefly,  it  can  be  argued  in  the  following  way.  Like  as 
in  a  two-dimensional  compressible  flow  the  advection  of 
floats  that  are  not  neutrally  buoyant,  gives  rise  to  clustering 
(Falkovich  et  al.,  2001).  The  floats  tend  to  dispose  along 
an  ocean  front  or  to  reach  a  convergence  zone.  The  float 
clustering  evolution  strongly  depends  on  the  structure  of  ver¬ 
tical  velocity  which  is  not  in  general  reproduced  accurately 
by  numerical  ocean  models.  Therefore  a  drift  and  spread  of 
a  cluster  composed  from  the  real  floats  differ  considerably 
from  those  of  a  cluster  of  the  synthetic  particles  even  if  the 
non-divergent  velocity  is  modeled  quite  accurately. 

Thus,  we  shouldn’t  expect  the  long-term  Lagrangian  pre¬ 
dictability  for  a  high-resolution  ocean  model  because  the  real 
drifters  and  synthetic  particles  have  different  asymptotic  be¬ 
havior. 
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2.2  Model  prediction  error 

Let  P  quasi-Lagrangian  floats  drogued  at  the  vertical  horizon 
z  display  the  spatio-temporal  variability  of  circulation  in  an 
area  of  interest.  The  model  may  reproduce  trajectories  of 
synthetic  floats  (particles)  to  be  compared  to  the  real  drifters. 

Let  the  position  and  velocity  of  some  drifter  and  synthetic 
particle  used  in  the  comparison  be  represented  by 

Z(t)  =  [X(t),  Y(t),  U(t),  V(t)}’  ,  (1) 

and 

Zpart(r)  =  {Xpart(0,  Tpart(r),  Cpart(r),  Vpart(r)}'  ,  (2) 

respectively,  the  superscript  V  denotes  the  transpose.  Then, 
the  model  prediction  error  be  represented  by  a  4-dimensional 
vector 

SZ(t)  =  {8X(t),  SY(t ),  8 U(t),  <5V(r)}?  ,  (3) 

with 

8X(t)  =  X(t)  -  Xpart(f),  8Y(t)  =  Y(t)  -  ypart(f), 

SU(t)  =  U(t)  -  Lpart(r),  8V(t)  =  V(t)  -  Lpart(t)  . 

Because  the  positioning  drifter  buoys  by  navigation  sys¬ 
tems  is  quite  accurate  with  error  less  than  200  m,  we  take 
that  differences  in  initial  positions  of  the  drifter  and  particle 
equal  to  zero  at  to-  That  is  not  true  for  difference  between 
their  velocities.  Thus  <5Z(fo)={0,  0,  8Uq,  SVo}?. 

Since  the  MPE  is  a-priori  stochastic,  the  Lagrangian  pre¬ 
dictability  is  effectively  determined  using  the  multi-variable 
PDF,  P(8X ,  8Y,  8U,  8  V.  t )  and  the  statistical  moments 

Lq={\8Z\«),  q  =  l,...,Q.  (4) 

The  MPE 

J  =  [\\8Z,A8Z\\a)  ,  (5) 

1 1  ...  1 1 4  is  the  4-dimensional  Eucledian  norm,  A  is  the  weight 
matrix,  determines  a  prediction  time  scale  r  from  the  condi¬ 
tion  that 

J(r)  <  e2  ,  (6) 

as  a  time  interval  within  which  this  RMS  error  is  less  than 
the  accepted  prediction  accuracy  e  (the  tolerance). 

In  general  the  use  of  Eqs.  (4)  and  (6)  sometimes  induces 
problems  in  physical  interpretation  of  the  MPE  growth  since 
different  dynamical  regimes  of  the  error  may  correspond  to 
the  same  power  decay  law  for  Lq  (Rangarajan  and  Ding, 
2000).  In  other  words,  it  is  difficult  (or  even  impossible) 
to  distinguish  different  physical  mechanisms  inducing  the 
anomalous  dynamical  regimes  through  Lq  in  practical  ap¬ 
plications. 


2.3  Irreversible-skill  time 

The  IT  (fin-ey)  has  been  proposed  for  measuring  the 
prediction-skill  without  the  Gaussian  distribution  assumption 
for  MPE  (Ivanov  et  al.,  1994;  Ivanov  and  Margolina,  1999) 
on  the  base  of  the  first-passage  time  (Gardiner,  1985).  This 
prediction  time  scale  is  determined  as  the  time  period  when 
J  exceeds  e2  for  the  first  time.  The  returning  prediction-skill 
that  is  referred  to  J  is  smaller  than  e2  after  IT  is  excluded. 

In  comparing  with  the  MPE,  two  different  statistics  for 
IT  can  be  defined  with  the  initial  error  8Z o  and  toler¬ 
ance  level  e:  (a)  PDF  (probability  density  function)  of  IT, 
W(r,  <5Zo,  e),  and  (b)  PSP  (probability  of  successful  pre¬ 
dictions),  (F(to.  8Zo,  e,  t—to)),  which  is  the  probability  that 
Tjrrev  is  larger  than  the  given  time  period  (t—to). 

The  PSP  satisfies  the  so-called  Pontryagin-Kolmogorov 
equation  (Chu  et  al.,  2002a,  b,  c).  The  statistics  are  con¬ 
nected  by 

t—to 

F(to,  <5Zo,  e,  t  —  to)  =  1  —  J  W(r,  SZq,  e)  dr  .  (7) 

o 

PSP  is  useful  to  verify  the  prediction-skill  in  practical  ap¬ 
plications  (Ivanov  et  al.,  1994;  Ivanov  and  Margolina,  1999) 
because  it  is  more  computational  feasible  than  the  PDF  of  IT. 

Additionally,  the  statistical  moments  of  IT  also  estimate 
the  practical  model  prediction  skill.  They  are  easily  calcu¬ 
lated  by 

OO 

Tk(SZo,  e)  =  k  J  F(to.  8Zq,  e,  t  —  to)(t  —  to)k~idt  .  (8) 

to 

The  mean  and  variance  of  IT  are  calculated  as 

("^irrev)  =  ^1  >  (9) 

(^rrev)  =  2r2  ~  *1  ■  (10) 

The  use  of  IT  allows  analytical  estimating  the  prediciton  skill 
through  the  statistical  moments  which  are  a  measure  of  any 
deviations  of  PDF  from  the  Gaussian  shape  (Ivanov  et  al., 
1994;  Chu  et  al.,  2002a,  b,  c).  Besides,  robust  algorithms 
for  the  estiamtion  of  PDF  and  PSP  from  numerical  modeling 
are  developed  by  Ivanov  et  al.  (1999)  and  Chu  et  al.  (2002a, 
2003c)  even  when  the  prediciton  skill  is  controlled  by  several 
variables  of  model  input. 

2.4  Singular  vector  approach 

This  approach,  commonly  used  to  estimate  large  atmospheric 
model  prediction-skill,  is  valid  only  for  small  MPEs  (Palmer, 
2000).  The  MPE  dynamics  is  the  linearization  with  respect 
to  the  reference  trajectory  Zref  described  by  the  following 
tangent-linear  model  (Lorenz,  1965;  Lacarra  and  Talagrand, 
1988  and  others)  within  the  reference  period  (T—to). 

~j~8Z(t)  —  Ji  [Zref (/)]  8Z  , 
at 

8Z(t0)  =  O  [Zreffo)]  tC[tQ,T],  (11) 
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where 

,  ,  3W/(Zref) 

is  the  Jacoby  matrix;  H  is  the  dynamical  operator  of  a  pre¬ 
diction  model. 

The  formal  solution  of  tangent-linear  Eq.  (1 1)  is  given  by 

SZ(t )  =  G(t,  t0)SZ(t0)  , 

where  G(t,  to)  is  the  transition  operator. 

The  eigenvalues  and  eigenfunctions  of  the  operator 

L  =  G'G  , 


are  called  the  singular  values  and  vectors  (“the  optimal  per¬ 
turbations”),  respectively.  The  directions,  along  which  the 
growth  of  SZ(to)  to  <5Z(r)  is  most  fast,  are  identified  by  the 
singular  vectors.  This  growth  leads  to  the  maximization  of 
the  amplification  factor 

(G-SZ,  G  SZ )  ( L  SZ ,  S Z) 

F  =  - - =  - - - - -  .  (12) 

(SZ,  SZ)  (SZ,  SZ) 

Geometrically  the  equation 


(L  ■  SZ,  SZ)  =  1  , 


(13) 


defines  a  hyper-ellipsoid  (ellipsoid  of  prediction  uncertainty) 
with  the  smallest  axis  maximizing  the  amplification  factor  F . 
This  axis  is  associated  with  the  largest  singular  value.  The 
fastest  growing  perturbations  are  the  singular  vectors  corre¬ 
sponding  to  the  largest  singular  values.  Knowledge  of  the 
leading  singular  vectors  allows  computating  the  amplifica¬ 
tion  factor  and  e-folding  predictability  time  (Palmer,  2000). 

The  singular  vectors  and  values  have  not  been  yet  ap¬ 
plied  as  a  measure  of  Lagrangian  model  predictability,  be¬ 
fore.  Thus,  two  problems  will  be  in  focus  of  the  present 
study.  First,  we  would  like  to  develop  a  simple  approach 
to  compute  the  singular  values  and  vectors.  Second,  we  need 
to  estimate  the  reference  period  from  the  real  oceanographic 
observations. 


2.5  Ensemble  realizations 


It  was  pointed  out  in  numerous  studies  (for  further  discus¬ 
sions  see  Palmer,  2000)  that  error  growth  may,  in  fact,  be 
highly  dependent  on  the  local  properties  along  the  attractor. 
Thus  our  goal  is  to  determine  “the  global”  mean  prediction 
skill  (where  “global”  average  is  to  be  understood  as  a  mean 
over  the  circulation  attractor).  Another  reason  to  average 
along  the  attractor  is  because  in  oceanographic  practice  the 
numbers  of  drifters  and  comparative  events  when  the  drifter 
and  particle  positions  are  compared,  are  not  large.  For  exam¬ 
ple,  only  55  drifters  moved  in  the  Gulf  of  Mexico  within  the 
time  period  under  investigation.  Therefore,  we  need  to  con¬ 
struct  an  ensemble  for  statistical  significant  estimations  from 
such  poor  observations. 

Ideally,  the  averaging  may  be  constructed  in  the  following 
way.  Utilizing  observations  from  P  Lagrangian  drifters  we 


firstly  reconstruct  the  spatial  structure  of  circulation  in  the 
area  of  interest  and  compare  it  with  computed  by  a  model  in 
a  prediction  metric.  For  simplicity  let  us  take  the  traditional 
RMS  error  written  as 

llsz2(t,  Zref(tQ))\\  = 

\\  II theor 

J  SZ2 (t)fi(SZ,  Z ref (t0),  t)dSZ  ,  (14) 

where  f\  is  the  PDF  for  a  model  prediction  error  over  the 
area  of  interest.  We  specially  stressed  in  Eq.  (14)  that  the 
model  prediction  skill  in  general  must  depend  on  initial  po¬ 
sition  Zref(fo)  along  the  attractor  trajectory. 

The  mean  error  is  determined  through  averaging  Eq.  (14) 
over  the  attractor  trajectory 


(15) 


[t. 


Z  (to)} 


I  B  [Zref(fo)]  h  [ZrefUo)]  dZre f 

theor 


where  /2  describes  the  distribution  of  comparative  events 
along  the  attractor,  B  is  the  weight  accounting  the  hetero¬ 
geneous  feature  of  the  attractor. 

In  practice  Eq.  (14)  is  transformed  to 

[[sz2(t)))  =  ^SZ2p[t,Z^(t0)]  (16) 

p= 1 

It  was  pointed  out  in  Vapnik  (1983)  and  Chu  et  al.  (2003a,  b) 
that  comparative  events  may  be  selected  so  that  for  the  em¬ 
pirical  average  (( SZ2(t )))  with  the  probability  larger  than  /, 
homogeneously  (without  excurses)  tends  to  the  theoretical 
one  if  P  — > oo,  i.e. 

Prob  sup  [[SZ2  [r,  Zref(ro)])) 

-(( SZ2(t )))  >x}->0  (17) 

This  explains  why  in  practice  impact  of  the  spatial  hetero¬ 
geneity  on  the  estimations  of  model  prediction  skill  may  be 
strongly  reduced  even  for  quite  heterogeneous  drifter  cover¬ 
ing  an  area  of  interest. 

The  most  dangerous  procedure  is  the  averaging  over  an 
attractor  trajectory.  In  the  present  study  we  will  apply  the 
following  procedure 

(sz2(r)}  =  J  /2(Zref)  ((<5Z2  [; t ,  Zref(fo)]))  dZKf  (18) 

where  B  —  1. 

Operator  (18)  represents  a  two-step  averaging  process:  (a) 
over  each  pair  of  drifter-particle  and  (b)  along  the  attractor 
trajectory.  In  the  present  study,  /2  is  assumed  to  be  the  con¬ 
tinuous  uniform  distribution  function,  i.e.  the  model  predic¬ 
tion  skill  is  stationary  or  at  least  very  slowly  variable  within 
the  reference  period. 

This  is  a  quite  strong  hypothesis  and  may  not  be  relaxed 
here,  because  the  number  of  drifter  buoys  equals  only  to  55. 
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3  Structure  of  model  attractor 


3.1  Model  circulation 


To  show  the  verification  of  the  Lagrangian  predictability  as 
an  example  the  Gulf  of  Mexico  real  time  nowcast/forecast 
system  is  taken.  This  system  (Kantha  et  al.,  1999;  Schaudt 
et  al.,  2001;  see  also  Kantha  and  Clayson,  2000)  is  based 
on  the  University  of  Colorado  version  of  the  Princeton 
Ocean  Model  (POM)  with  1/12°  horizontal  resolution  and 
24  er—  levels.  The  water  inflow  through  the  Yukatan  Channel 
has  been  prescribed,  the  outflow  through  the  Strait  of  Florida 
evolved  with  the  model  physics.  The  model  reaches  quasi 
steady  regimes  after  10  year  integration.  High-resolution 
bathymetry  (Fig.  la),  the  real-time  altimetric  sea  surface 
height  (SSH)  anomalies  derived  from  NASA/CNES  TOPEX 
and  ESA  ERS-1/2  altimeters,  and  composite  MCSST  data 
derived  from  NOAA  AVHRR  assimilated  into  the  model  in 
a  continuous  data  assimilation  mode  are  used  to  produce 
a  nowcast  and  a  four-week  forecast.  Kantha  and  Clayson 
(2000)  found  that  forecast  retains  considerable  skill  to  about 
1-2  weeks,  beyond  which  the  forecast  begins  to  deviate  from 
reality. 

Since  the  Lagrangian  drifters  are  drogued  at  50  m  depth, 
the  modeled  horizontal  velocity  field  at  the  same  depth  is 
used  for  the  model-data  comparison  in  1998,  throughout  this 
paper.  The  model  successfully  reproduces  the  typical  circu¬ 
lation  patterns  (the  Loop  Current  and  eddy  shedding)  in  the 
Gulf  of  Mexico  (Fig.  lb,  c).  The  anti-cyclonic  eddies  are  reg¬ 
ularly  pinched  off  from  the  Loop  Current,  which  enters  the 
eastern  Gulf  of  Mexico  through  the  Yucatan  Channel  and  ex¬ 
its  through  the  Strait  of  Florida.  The  Loop  Current  meanders 
and  eddies  (most  distinguished  mesoscale  features)  are  pre¬ 
dicted  with  radii  less  than  100  km  and  velocities  larger  than 
100  cm/s. 


3.2  Toroidal  and  poloidal  EOF  decomposition 


The  model  circulation  at  50  m  depth  was  analyzed  as  a  2- 
D  divergent  flow.  Helmholtz  (Morse  and  Feshbach,  1953) 
or  Moffat-Zeldovich  (Moffat,  1978;  Zeldovich  et  al.,  1985) 
decomposition  may  be  applied  in  this  case.  Our  study  uses 
the  MZD  to  construct  appropriate  EOFs  for  the  analysis  of 
the  spatio-temporal  structure  of  the  model  attractor  at  50  m 
depth.  Other  oceanographic  applications  of  MZD  can  be 
found  in  Eremeev  et  al.  (1992),  Lipphardt  et  al.  (2000)  and 
Chu  et  al.  (2003a,  b)  selected  as  examples. 

Let  us  briefly  explain  the  MZD.  It  represents  3-D  velocity 
\u_,  ?/;)  at  any  vertical  horizon  z=const  as 


( t ) 

!/J_  =  U 

„(C  _ 


,(P) 


=  V  x  (£'!') . 
w  =  —V  ■  V<t> , 


u(P)  =  V  x  V  x  (*<$>) , 


(19) 

(20) 
(21) 


where  k  is  the  unit  vector  perpendicular  to  the  surface 
z=const;  uj_=(u,  u)  and  w  are  the  horizontal  and  vertical 
velocity  components,  respectively;  V=(VA,  V  v )  is  the  gra¬ 
dient  operator.  The  T  (toroidal)  and  <J>  (poloidal)  potentials 


Fig.  1.  High-resolution  bathymetry  utilized  by  the  model  (a)  and 
circulation  patterns  in  the  Gulf  of  Mexico  for  two  different  dynam¬ 
ical  regimes  in  1998  on  (b)  Julian  Day  260,  and  (c)  Julian  Day  320. 
The  bold  line  in  (a)  is  50  m  isobath.  Transition  from  one  regime  to 
another  is  characterized  by  energy  cascade  from  low-  to  high-order 
EOFs. 
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generate  the  toroidal  ( u ' 1 1 )  and  poloidal  (u'p>)  velocities, 
respectively. 

Knowledge  of  the  potentials  rp  and  <J>  allows  computing 
EOFs  separately  for  each  of  these  potentials:  T,,  (x,  y,  z), 
7i=l , . . . ,  N  (toroidal  modes)  and  <t>„,  (x,  y,  z),  ;«=1 , . . . ,  M 
(poloidal  modes). 

The  following  scheme  was  chosen  for  practical  compu¬ 
tations  of  the  EOFs.  Let  T— to  be  the  reference  period 
during  which  the  model-data  comparison  takes  place.  In 
our  case  this  time  period  can  be  chosen  as  not  larger  than 
1  year.  Then,  the  horizontal  velocity  u±=(u(T\  it' />l  ) 
and  the  potentials  ('P,  <f>)  were  represented  as  combina¬ 
tions  of  stationary  within  the  reference  period  currents 
(SPCs)  u±=  (u(T\  u{P)) ,  Tp  <i>  and  pulsation  components 

u\=(u,(T\  i i,<p>y  vp',0': 


4r(x,  y,  z,  t)  —  ^(x,  y,  z)  +  ^'(x,  y,  z,  t), 
u(p  =  «(P)  +u’{P\ 

d>(x,  y,  z,  t )  =  0(x,  y,  z)  +  ^'(x,  y,  z,  t), 
u{T)  =  u(T)+u,(T). 


(22) 

(23) 


Two  sets  of  EOFs  which  are  the  eigenfunctions  of  a 
large  18  200x18  200  matrix,  were  computed  for  P'  and  <P' 
through  the  algorithm  of  Penenko  and  Protasov  (1978).  To 
compute  EOFs  we  used  the  same  grid  as  in  the  hydro- 
dynamic  model  and  circulation  velocity  through  each  two 
days  of  model  calculations.  Numerical  experiments  demon¬ 
strated  that  for  the  reference  period  larger  than  8  months  the 
computed  EOFs  were  practically  not  sensitive  to  increasing 
length  of  this  period. 

The  spatial  inhomogeneity  and  temporal  intermittency  of 
the  circulation  attractor  are  analyzed  through  the  spatial 
structure  of  EOFs  and  behavior  of  mode  amplitudes  A„ 
and  Bm  calculated  by  the  inner  products  T/i  =  (P'  ■  P„)  and 
Z?,„  =  (< f>'  •  <P,„),  respectively. 


3.3  Spatial  structure  of  model  attractor 


3.3.1  Stationary  part 


Both  the  toroidal  (ii(r))  (Fig.  2a)  and  poloidal  (ii(P)) 
(Fig.  2b)  SPCs  have  simple  large-scale  topology.  The 
toroidal  circulation  explicitly  dominates  over  the  poloidal 
one.  This  indicates  the  quasi-geostrophic  nature  of  SPC  pat¬ 
terns.  A  measure  of  toroidal  circulation  impact  is  the  com¬ 
pressibility  factor  (Schumacher  and  Eckartdt,  2002) 


(((Vm_l)2» 


(24) 


i.e.  the  ratio  of  the  mean  square  divergence  versus  the  mean 
square  velocity  gradient.  Here,  the  double  brackets,  also  as  in 
Eq.  (14),  denote  averaging  over  an  area  of  interest.  The  com¬ 
pressibility  factor  of  the  model  SPC  is  about  0.05.  However, 
the  poloidal  SPC  (m(/>))  (Fig.  2b)  clearly  shows  the  existence 
of  the  convergence  and  divergence  zones  in  the  circulation 
patterns. 


Lagrangian  predictability  of  high-resolution  regional  models 


Fig.  2.  Structures  of  the  toroidal  (a)  and  poloidal  (b)  SPCs. 


3.3.2  Pulsation  component  of  circulation 

The  spatial  structures  of  the  non-stationary  part  are 
represented  by  various  toroidal  (n=l,...,N)  and 
poloidal  velocity  modes  (snapshots) 

{hJ^Vx [£47„]}  and  {ul^—V  ■  k<i>m}.  The  most  striking 

(T) 

features  of  {un  }  are  (a)  the  multi-eddy  structures  with 
smaller  eddy  scales  for  high-order  modes,  (b)  the  Loop 
Current  occurring  only  in  the  first  toroidal  mode  u\  but 
not  in  other  modes  and  (c)  small-scale  eddies  contributing  in 
the  low-order  modes  as  well.  For  example,  the  first  toroidal 
velocity  mode  u\  represents  a  combination  of  the  Loop 

Current  and  multi-eddy  structures  (Fig.  3a).  The  second 

(T) 

toroidal  velocity  mode  u2  (Fig.  3b)  shows  the  existence 
of  multi-eddy  structures  such  as  dipoles.  The  most  striking 
feature  of  {«,„  }  is  that  the  convergence  and  divergence 
zones  at  the  mesoscales  in  the  deep  basin  are  represented 
only  by  the  poloidal  velocity  modes  with  order  higher  than 
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Fig.  3.  Low-order  toroidal  [n=l  (a),  n=2  (b)]  and  poloidal  [m=l  (c),  m=2  (d)]  EOFs. 


two.  The  first  (n\F))  (Fig.  3c)  and  the  second  (Fig.  3d) 
poloidal  velocity  modes  display  the  large-scale  motions  as 
well  as  the  upwelling  and  downwelling  in  the  coastal  zone 
of  the  Gulf  of  Mexico. 

The  fractions  of  energy  for  each  type  of  EOFs  grow  rapidly 
with  the  increase  of  the  mode  numbers  (N,  M)  (Fig.  4a,  b). 
Here,  more  than  90%  of  the  toroidal  and  99%  of  the  poloidal 
kinetic  energy  lie  in  the  first  ten  toroidal  and  two  poloidal 
modes,  respectively.  The  first  toroidal  mode  contains  more 
than  45%  of  the  total  energy  of  toroidal  motions. 

3.4  Temporal  structure  of  model  attractor 

Temporal  variations  of  the  first  ten  toroidal  mode  amplitudes 
{Ai(0,...A  io(f)}  (Fig.  5)  show  various  degrees  of  com¬ 
plexity,  which  increases  in  the  higher-order  modes.  Fairly 
large  excursions  from  the  mean  are  also  noted.  The  temporal 
variability  scale  of  the  second  (Fig.  5b)  and  third  (Fig.  5c) 


modes  exceeds  the  temporal  variability  scale  for  the  first 
mode  (Fig.  5a).  Amplitudes  of  the  sixth  (Fig.  5f),  ninth 
(Fig.  5i),  and  tenth  (Fig.  5j)  modes  demonstrate  irregular 
behaviors.  Temporal  variations  of  the  first  two  poloidal- 
mode  amplitudes  {B\  ( t ),  /T(f)}  show  the  similar  complexity 
(Figs.  5k  and  51). 

The  pulsations  of  mean  toroidal  {{E[0 r))  and  poloidal 
(j\Epo\j)  energies  also  demonstrate  complex  irregular  (abrupt) 
behavior  (Fig.  6a,  b).  The  ((£( or))  varies  between  240  and 
600cm2/s2  with  energy  peeks  occurring  at  Julian  Days-150, 
-210,  -300,  and  -340  of  1998  (Fig.  6a).The  relative  variations 
of  are  much  stronger,  between  10  and  80cm2/s2, 

with  explicit  energy  peeks  at  Julian  Day-230,  -245  and  -360 
of  1998  (Fig.  6b).  The  behavior  of  [(E'tOT))  is  de-locked  with 

))■ 


the  variations  of  ( ( EL 
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Number  of  Modes 


Fig.  4.  The  fractions  of  toroidal  (a)  and  poloidal  (b)  energies.  Ar¬ 
rows  indicate  the  fraction  of  energy  corresponding  to  1,  10  and  20 
EOF  decomposition. 


The  complexity  in  the  circulation  attractor  behavior  can  be 
also  analyzed  using  the  four  non-dimensional  ratios: 


_  Ogroi))  „  \\  po1// 

m~  «£tor»  ’  1l2-((E'tm))- 


-  -  -  - (M  onrl 

'n~  «Gor»  ’  ''-_{{£;or}>’  ^4_((£pol»’ 

global  entropy  (Aubry  et  al.,  1991).  Here,  the  mean 
total  ((E)),  toroidal  ((£tor)}  and  poloidal  ((£p0i))  ki¬ 
netic  energies  are  computed  as  ((E))  —  ((Etm))  +  ((_Epoi)), 


netic  energies  are  computed  as  ((E))  —  ((Etor))  +  ((_Epoi)), 
((£tor»=ptor))  +  ((£t'or))  and  [[Epol))  =  [[Epol))  +  ({E'pol)), 
respectively. 

The  ratios  ?;i  and  772  vary  between  2%  and  12%  (Fig.  7a) 
and  between  2%  and  21%  (Fig.  7b),  respectively.  The  be¬ 
havior  of  both  the  ratios  are  locked.  That  indicates  the  same 


forcing  mechanism  for  the  Gulf  of  Mexico  circulation  at  the 
large  and  meso-scales. 

The  ratios  773  and  774  fluctuate  from  42%  to  73%  (Fig.  7c) 
and  from  3%  to  27%  (Fig.  7d),  respectively.  We  see  that  the 
pulsations  of  the  poloidal  component  contribute  very  strong 
into  the  total  poloidal  circulation.  That  stresses  our  conclu¬ 
sion  about  impossibility  of  the  global  Lagrangian  predictabil¬ 
ity  for  high-resolution  regional  models. 


200  250  300  350  200  250  300  350 

Model  Days  Model  Days 


Fig.  5.  Temporal  evolution  of  mode  amplitudes  A\, . . . ,  Ap)  at 
(a) . (j)  and  B{ ,  B2  at  (k),  (1). 
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Fig.  7.  Temporal  evolution  of  non-dimensional  ratios  (a),  r] 2  (b), 
m  (c)  and  tj 4  (d). 


Two  global  entropies,  toroidal  (Stor)  and  poloidal  i.Spoi) 
can  be  calculated  for  the  toroidal  and  poloidal  movements, 
respectively,  as 


N 

Star  =  -£>„  log  Pn/  log  N  , 
n=  1 


M 

Spol  =  -  Pm  lo8  Pm/k>&  M  > 


m=  1 


N 

Pn  =  A2n/J2Al’  ^ 

n— 1 
M 

P*m  =  4/  £  Bm  ■  C26) 

m=  1 


The  entropy  .S't()l-  or  .S'po]  is  0  (total  order)  if  and  only  if  the  first 
EOF  mode  dominates.  If  the  energy  is  distributed  among  the 
modes  equally  (total  disorder),  the  entropy  takes  its  maxi¬ 
mum  value,  namely  1 .  Thus,  the  global  entropies  are  a  mea¬ 
sure  of  degree  of  order  or  disorder  in  behavior  of  a  model 
attractor. 

The  global  toroidal  and  poloidal  entropies  with  mode  trun¬ 
cation,  N—M— 16,  and  N — M  —50,  are  illustrated  in  Fig.  8. 
Both  these  entropies  display  intermittent  characteristics  by 


Fig.  8.  Temporal  behavior  of  the  toroidal  (a,  b)  and  poloidal  (c,  d) 
entropies  computed  for  16  (a,  c)  and  50  (b,  d)  EOFs. 


successive  redistribution  of  energy  initially  concentrated  in  a 
few  first  low-order  modes  (Stor=0.2  on  Julian  Day-260),  be¬ 
tween  other  low-  and  high-order  EOFs.  Such  intermittency 
is  more  evident  in  the  poloidal  field  (Fig.  8c,  d)  than  in  the 
toroidal  one  (Fig.  8a,  b).  As  a  result  displayed  in  Fig.  9a-d  is 
the  existence  of  irregular  (probably,  stochastic)  pathways  of 
divergence  and  convergence  zones. 

4  Lyapunov  dimension  of  model  attractor 

The  global  prediction  feature  of  the  model  attractor  can  be 
identified  through  different  ergodic  measures,  such  as  Lya¬ 
punov  dimension  ( /)lyap),  Kolmogorov-Sinai  entropy  and 
others  (Zaslavsky,  2002).  Exact  determination  of  them  for 
a  high-resolution  oceanographic  model  is  a  very  hard  task. 
However,  coarse  estimations  of  Z)Lyap  are  possible  through 
the  approach  originally  developed  by  Syrovich  (1989). 

The  minimum  mode  truncation  number  required  for  accu¬ 
rate  representing  the  circulation  as  measured  by  the  energy 
norm  is  treated  as  the  intrinsic  dimension  Dmt  of  phase  space 
generated  by  the  basis  { 4',, (x ,  y)}  and  {T>,„(x,  y)}.  This  di¬ 
mension  is  identified  as  the  number  of  modes  requested  so 
that  the  captured  energy  is  at  least  90%  of  the  total  one  and 
so  that  no  neglected  mode,  on  average,  contains  more  than 
1%  of  the  energy  lying  in  the  most  energy  mode  (Syrovich, 
1989). 

It  was  obtained  in  numerous  studies,  here  we  cite  only  Sy¬ 
rovich  (1989)  and  Aubry  et  al.  (1991),  that  the  intrinsic  and 
Lyapunov  dimensions  are  connected  through  the  following 
empirical  relation 

^DLyapj  =  Dint  -  1,  (27) 

where  [,DLyap]  is  the  integer  part  of  Z)Lyap. 
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Fig.  9.  Model  reproduced  poloidal  circulation  in  1998  on  Julian  Day  200  (a),  Julian  Day  235  (b),  Julian  Day  300  (c),  and  Julian  Day  360  (d). 


Since  the  circulation  is  decomposed  into  the  toroidal 
and  poloidal  parts,  two  different  intrinsic  dimensions  and 
naturally  two  Lyapunov  dimensions  should  be  introduced 
as  D^aP)  and  ^poiaP)  f°r  t'le  toroidal  and 

poloidal  velocity  components,  respectively.  The  analysis  of 
fractions  of  toroidal  and  poloidal  energies  (see  Fig.  4a,  b)  re¬ 
sults  in 


D" 


16; 


D 


15  <  /j'fip  <  16 


poi  “  2 ; 
n 

pol 


1  <  £>PI,ap  <  2 


(28) 

(29) 


Thus,  both  the  estimations  for  D^ap  and  £>pjap  uniquely  in¬ 
dicate  the  irregularity  of  behavior  of  the  model  attractor.  That 
can  lead  to  the  exponential  growth  of  the  initial  error  on  the 
attractor. 


5  Statictics  of  the  Lagrangian  predictability 


The  present  section  analyzes  the  predictability  in  general 
with  uncertainty  in  the  evolution  law  (the  governing  equa¬ 
tions)  and  partially  in  initial  conditions.  Several  physical 
factors  can  limit  the  model  prediction  in  the  Gulf  of  Mex¬ 
ico.  First,  the  water  inflow  through  the  Yukatan  Channel 
and  the  outflow  through  the  Strait  of  Florida  are  described 
only  approximately.  The  prediction  error  induced  by  inac¬ 
curacy  in  open  boundary  conditions  can  considerably  reduce 
the  model  prediction  time  (Chu,  1999;  Jiang  and  Malanotte- 
Rizzoli,  1999).  Second,  other  sources  of  prediction  errors 
are  unresolved  motions  which  are  either  parameterized  in  the 
model  or  neglected.  Comparison  between  the  drifter  trajecto¬ 
ries  and  modeled  synthetic  particle  trajectories  (Toner  at  al., 
2001)  demonstrates  that  after  filtrating  the  effect  of  30-hour 
tides,  individual  drifter  trajectories  are  reproduced  by  the 
model  with  reasonable  accuracy  within  10  days.  The  present 
study  utilizes  the  same  drifter  observations  and  numerical 
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Fig.  10.  Spaghetti  of  55  drifting  buoys  used  in  the  analysis. 


model.  However,  in  comparing  with  Toner  et  al.  (2001)  and 
Kuznetsov  et  al.  (2002)  the  primary  focus  of  our  study  is  on 
statistical  laws  of  the  model  prediction  error  growth. 

5.1  Drifter  observations 

Fifty  five  satellite-tracked  sonobuoys,  manufactured  by  the 
Applied  Technology  Associates  and  deployed  by  the  Hori¬ 
zon  Marine,  have  provided  the  Lagrangian  data  used  for 
the  study.  A  nylon  drogue  at  50  m  depth  is  tethered  to  the 
buoy  hull  that  is  0.64  m  in  the  water  and  0.33  m  in  the  air. 
Louisiana  State  University’s  Coastal  Studies  Institute/Earth 
Scan  Laboratory  tracked  the  buoys  for  the  Horizon  Marine. 
The  drifter  trajectories  used  in  the  present  study  are  given  in 
Fig.  10.  The  particle  trajectories  are  obtained  through  time- 
integration  of  the  modeled  horizontal  velocity  field  at  50  m 
depth  with  the  bi-cubic  spline  interpolation  in  space  and  the 
second-order  interpolation  in  time. 

5.2  Model  prediction  error 

The  differences  of  position  (8X,  8Y)  and  velocity  (8U .  8  V ) 
for  any  drifter-synthetic  particle  pair  are  stochastic.  The  en¬ 
semble  average  of  the  velocity  difference  over  the  trajecto¬ 
ries  (( SU )  ,  (SV))  shows  weak  model  bias.  For  example. 


the  bias  (SU)  has  a  minimum  value  of  -2.5  cm/s  on  Julian 
Day-4  in  1998  and  a  maximum  value  of  5  cm/s  on  Julian 
Day-18  in  1998  (Fig.  11a).  However,  the  absolute  RMS  er¬ 
rors  of  velocity  components  (over  the  trajectories)  are  not 
small.  For  example,  ( SU2)1 ^  fluctuates  between  34.5  cm/s 

and  44.4  cm/s,  and  (<5  V2} 1  ^  fluctuates  between  3 1 .7  cm/s  and 
38.0 cm/s  (Fig.  lib).  The  variances  of  prediction  error  for 
both  the  drifter  velocity  components  become  practically  con¬ 
stant  after  4  days  from  the  start  of  particle  movement.  In  gen¬ 
eral  large  stationary  variances  of  MPE  indicate  that  the  error 
evolutes  nonlinearly  (Nicolis,  1992;  Chu  et  al.,  2002a). 

Let  us  discuss  how  the  position  model  error  characterizing 
differences  between  the  drifter  and  particle  positions  evolves. 
We  will  determine  its  statistics  by  the  first  ten  moments 

Lq(t)  =  ^SZ9]j  4  =  1,...,  10,  (30) 

where  8Z={8X,  8Y}' . 

It  was  pointed  out  by  Chu  et  al.  (2002a)  that  Eq.  (30)  holds 
a  power  law  for  q—  1  and  q—2.  Here,  we  check  the  existence 
of  a  more  universal  power  law  for  4=1,  . . . ,  10 

Lq(t)  ~  tYq  (31) 

and  specify  the  power  exponent  yq. 
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Fig.  11.  Mean  (a)  and  variance  (b)  of  difference  between  buoy  and  particle  velocities. 


Fig.  12.  The  first  ten  moments  of  MPE.  The  solid  lines  represent 
the  power  laws  which  approximate  computation  results  displayed 
by  dots,  crossings,  triangles  and  squares. 


If  Yq=ml  2,  (32) 

the  PDF  of  MPE  has  a  self-similar  form.  MPE  described  by 
relation  (31)  with  (32)  is  a  Gaussian  (non-Gaussian)  process 
if  \x—  1  (/xyM)  (del  Castillo-Negrete,  1998).  Any  deviation 
from  relation  (32)  suggests  that  the  MPE  is  an  intermittent 
process. 


Fig.  13.  The  scaling  exponent  yq  as  measured  within  the  range 
12h<f<20  days.  Dashed  line  is  the  3g/4  law.  Solid  line  is  the 
best  fit.  Difference  between  the  3^/4  law  and  the  best  fit  observed 
for  the  moments  numbered  higher  than  4  is  explained  by  reducing 
accuracy  of  computational  results  from  finite  sampling  time  effects. 


Our  computations  display  a  power  dependence  of  Lq  (f )  on 
the  time  t  (12  h  <t< 20  days)  from  q— 1  to  t/=10  (Fig.  12): 

Lg(t)~tM'2  (33) 

with 

H  1.49  ±0.03, 
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Fig.  14.  Psx  =  Pi(SX,  0,  t)  (the  black 
triangles)  and  Pgy=Pi  (0,  ST,  t)  (the 
white  dots)  through  3h  (a),  6h  (b), 
12  h  (c),  24  h  (d),  48  h  (e)  and  96  h  (f). 
The  solid  lines  are  Gaussian  approxi¬ 
mations. 


8U,  SV  (m/s) 


(b) 
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Fig.  15.  PSU  =  P2(SU,0,t)  (the  black 
triangles)  and  Psv=Pl(0,  SV,  t)  (the 
white  dots)  through  3h  (a),  6h  (b), 
12  h  (c),  24  h  (d),  48  h  (e)  and  96  h  (f). 
The  solid  lines  are  Gaussian  approxi¬ 
mations. 


To  assess  the  possible  impact  of  finite  sample  size  effects 
on  estimating  /x,  we  apply  a  simple  test  that  was  proposed 
earlier  by  Tennekes  and  Wijngaard  (1972).  The  maximum 
error  for  estimating  n  is  about  0.02-0.03  for  the  high-order 
moments.  For  the  first  four  moments  Lq  (q  —  \ .  . . . ,  4)  the 
scaling  exponent  is  very  close  to  3g/4  (Fig.  13). 

The  power  growth  of  MPE  is  interpreted  in  terms  of 
the  normal  and  anomalous  diffusion  (del  Castillo-Negrete, 
1998).  The  cases  of  /x=  1  and  /x>  1  may  be  called  as  “the 
normal  diffusion”  and  “the  super-diffusion”  regimes  of  MPE 
growth,  respectively. 


5.3  PDF  of  MPE 

To  understand  how  the  statistics  of  SZ  departures  from  Gaus¬ 
sian  distribution,  the  following  shorted  PDFs  are  applied: 
Pi(SX,  8Y,  t)  and  P2(&U ,  SV ,  t) .  The  appropriate  cross- 
sections  of  these  PDFs  are  plotted  in  Figs.  14  and  15  display¬ 
ing  their  explicit  symmetry  and  isotropy.  Pi(SX,  ST,  1)  has 
a  non-Gaussian  shape  for  various  times  from  3  h  to  10  days. 
The  departures  from  Gaussian  behavior  and  the  presence  of 
tails  stretched  in  large  MPE,  are  both  evident  in  these  fig¬ 
ures.  At  times  longer  than  10  days  this  function  tends  to  a 
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Fig.  16.  Cluster  spread  along  the  trajectory  (a)  of  H17543  drifter.  Positions  1,  2,  3,  4,  5,  6  of  cluster  center  (a)  and  corresponding  to  them 
positions  of  the  cluster  (b),  (c),  (d).  (e),  (f),  (g)  are  given  through  every  8  h.  Initially  the  cluster  is  round  with  the  radius  equaled  to  1/12°  (b). 
The  solid  lines  are  the  cluster  boundaries  predicted  by  the  singular  vector  approach.  The  black  dots  are  particles. 


Gaussian  shape  because  of  correlations  for  a  pair  of  buoy- 
particle  decay.  The  shape  of  P2(8U ,  8V,  t)  practically  does 
not  change  with  time  (Fig.  15). 

Using  the  analogy  between  the  pair-particle  dispersion  and 
the  drifter-synthetic  particle  divergence,  it  is  reasonable  to 
assume  that  P\(8X ,  SY,  t)  holds  the  following  diffusion-like 
equation  (Boffeta  and  Celani,  2000): 


‘-m 


(34) 


where  SR= V 8X2  +  8Y2,  Kef{(R,  t )  ~  R 3/4  is  the 

Richardson-Obukhov  type  effective  diffusion  coefficient. 
The  solution  of  Eq.  (34)  is  written  as  (Klafer  et  al.,  1987) 


Psr  =  f-9|/2exp  (-C  ■  7?2//3/f)  , 


(35) 


here,  C  is  constant. 

Relation  (35)  correctly  describes  the  evolution  of  the 
shorted  PDF  of  MPE  for  the  small  and  intermediate  times. 
This  indicates  that  the  statistics  of  the  Lagrangain  prediction 
error  in  the  case  of  the  second  kind  of  predictability  (Lorenz, 
1984)  can  be  described  by  diffusion-like  equations.  Note  that 
the  model  error  is  scaled  within  the  time  intervals  from  12  h 
to  20  days.  These  time  scales  correspond  to  the  spatial  scales 
varying  from  20  km  to  300  km. 


5.4  Linear  tangent  model 

To  examine  applicability  of  linear  tangent  models  for  the 
analysis  of  Lagrangian  predictability  of  high-resolution  re¬ 
gional  models  we  study  the  dispersion  of  small-size  particle 
clusters  with  the  initial  round  shape  and  the  radius  equaled  to 
1/12°  along  the  drifter  trajectories  [Zref(f)=Z(t)]  and  com¬ 
pare  computation  results  with  predictions  from  the  singular 
vector  approach. 

If  a  linear  tangent  model  is  valid  within  the  reference  pe¬ 
riod  (T—to),  the  cluster  rapidly  gets  an  ellipsoidal  shape  con¬ 
served  within  this  period.  In  opposite  case,  the  cluster  is 
quickly  deformed  to  a  pancake-like  shape.  The  evolution  of 
the  cluster  shape  can  be  studied  in  two  ways:  through  direct 
modeling  and  by  the  singular  vector  analysis.  This  allows 
estimating  the  reference  period  for  the  real  circulation  model 
from  the  real  oceanographic  observations. 

The  results  discussed  above  say  that  the  Lagrangian  model 
skill  is  non-Gaussian  even  for  small  times.  Thus,  we  can  not 
expect  large  reference  periods.  However,  their  exact  estima¬ 
tions  are  of  considerable  interest  for  the  practical  predictabil¬ 
ity. 

In  general  we  found  that  linear  tangent  models  can  be  ap¬ 
plied  to  approximate  the  MPE  growth  within  very  short  time 
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Fig.  17.  Cluster  spread  along  the  trajectory  (a)  of  H17547  drifter.  Positions  1,  2,  3,  4,  5,  6  of  cluster  center  (a)  and  corresponding  to  them 
positions  of  the  cluster  (b),  (c),  (d),  (e),  (f),  (g)  are  given  through  every  15  h.  Initially  the  cluster  is  round  with  the  radius  equaled  to  1/12°  (b). 
The  solid  lines  are  the  cluster  boundaries  predicted  by  the  singular  vector  approach.  The  black  dots  are  particles. 


intervals  limited  by  several  days.  As  an  example,  two  drifter 
(H17543  and  H17547)  trajectories  are  selected  to  show  this 
(Figs.  16  and  17).  The  reference  periods  for  H17543  and 
H17547  drifters  were  estimated  as  8  and  30  h,  respectively. 
The  mean  reference  period  {(T—to))  was  equaled  to  21  h, 
i.e.  after  this  time  a  linear  tangent  model  breaks  down.  Thus, 
21  h  is  the  practical  limit  to  apply  a  linear  tangent  model  of 
MPE  in  the  Gulf  of  Mexico. 

The  analysis  of  evolution  of  all  clusters  used  to  estimate 
the  reference  period  also  demonstrate  that  although  a  clus¬ 
ter  is  quickly  deformed  to  a  pancake-like  shape,  the  cluster 
center  moves  along  the  drifter  trajectory  during  considerably 
larger  time  than  the  reference  period,  i.e.  the  local  model  drift 
was  a  quite  small. 

Because  the  mean  reference  period  should  be  limited  only 
by  21  h,  we  will  not  apply  linear  tangent  models  for  the  anal¬ 
ysis  of  Lagrangian  predictability  any  more. 

5.5  IT  statistics 

The  PSP  (see  Eq.  7)  with  zero  initial  errors,  F(t,  e),  is  calcu¬ 
lated  from  the  drifter  and  particle  trajectories  with  the  toler¬ 
ance  level  ranging  from  0.05°  to  1.25°.  F(t.  e)  for  four  dif¬ 


ferent  tolerance  levels  (0.25°,  0.5°,  0.75°,  and  1.25°)  clearly 
shows  non-Gaussian  distribution  (Fig.  18).  The  long  tails 
stretched  into  large  ITs  demonstrate  the  existence  of  the  long¬ 
term  correlations  between  the  drifters  and  synthetic  particles. 
We  see  that  individual  Lagrangian  trajectories  can  be  pre¬ 
dicted  within  20-25  day  period  even  if  the  mean  Lagrangian 
prediction  time  does  not  exceed  4-5  days.  Individual  predic¬ 
tions  with  abnormal  large  prediction  periods  were  called  the 
extremely  successful  predictions  (ESP)  (Chu  et  al.,  2002a). 

The  PSP  tails  (Fig.  18)  have  the  power  behavior 
for  the  long  IT  period  with  the  power  exponents  as 
(2.17±0.07,  1.98±0.09,  2.08±0.10,  1.77±0.25)  for  e  as 
(0.25°;  0.5°,  0.75°,  1.25°),  respectively.  That  identifies  the 
power  exponent  as  larger  than  2. 

If  the  power  exponent  is  larger  than  2,  the  mean  and  vari¬ 
ance  of  IT  exist  (Fig.  19).  Both  these  values  are  scaled  as 

(tirrev>~e“,  (^ev)  ~  ^  (36) 

with  aRG.21iO.ll  and  j6Ri2.17i0.15. 

That  indicates  non-Gaussian  feature  of  the  Lagrange  pre¬ 
diction  skill.  It  was  pointed  out  by  Popov  (1970)  that  a— 1 
and  /3= 2  for  a  Gaussian  statistics. 
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Fig.  18.  PSP  of  IT  calculated  for  different  levels  of  tolerance  e. 
Solid  and  dashed  lines  are  the  PSP  and  the  proposed  theoretical 
power  laws  to  approximate  it  for  large  times. 

The  mean  irreversible  time  can  also  be  easily  computed 
from  model  (34),  (35)  and  through  Popov’s  (1970)  method. 
The  theoretical  estimation  of  a  is  4/3.  We  see  that  the  best  fit 
of  the  power  exponent  estimated  as  1.21±0.11  agrees  with 
the  theoretical  value  very  well. 

6  Reconstruction  of  model  attractor 

It  was  demonstrated  above  that  mean  Lagrangian  prediction 
skill  of  high-resolution  oceanographic  models  is  not  very 
high,  up  to  4-5  days,  and  models  predict  the  real  drifter  tra¬ 
jectories  only  within  short  and  intermediate  time  periods.  A 
key  question  is:  may  the  spatial  structure  of  such  a  model 
attractor  be  reconstructed  from  the  rare  drifter  observations? 

An  approach  developed  in  Eremeev  et  al.  (1992)  and  Chu 
et  al.  (2003a,  b)  allows  reconstructing  velocity  field  from 
sparse  and  noisy  observations.  Accordingly  to  them  a  ve¬ 
locity  field  is  projected  onto  N  dimension  phase  space  gen¬ 
erated  by  some  basis  functions.  Then,  these  projections  (the 
reconstruction  coefficients)  are  estimated  from  P  observa¬ 
tions.  In  general  the  reconstruction  procedure  requests  that 
P>{\.5-2)N. 

We  will  apply  this  approach  for  estimating  the  Gulf  of 
Mexico  model  attractor  structure  from  15  drifter  buoys  if 
EOFs  computed  from  model  circulation  are  taken  as  the  basis 
functions. 

Although  the  EOF  decomposition  has  rapid  convergence 
by  the  energetic  norm,  the  accurate  approximation  of  spatial 
structure  of  circulation  needs  more  than  few  EOFs.  The  ac¬ 
curacy  of  the  EOF  truncation  is  represented  by  the  relative 
RMS  errors  of  the  toroidal  and  poloidal  potentials 

Xtor(L  N)  =  Ij^'lp1  H*'  — 4^11  (37) 

XPol(f,M)  =  Hep'll-1  | (38) 


where  the  Euclidian  norm  1 1  1 1  summarizes  over  all 

N 

points  of  the  model  computation  grid,  'T'v=  Y  A,,'!',,  and 

n= 1 
M 

tk 'M—  Y  Bm<t>n  are  the  pulsations  of  potentials  approxi- 

m=\ 

mated  by  N  toroidal  and  M  poloidal  modes,  respectively. 

Figures  20  and  21  show  the  temporally  varying  Xtor(L  N) 
and  Xpol  (t,  M)  with  six  different  mode  truncations  (1,  3,  6, 
10,  30,  50),  respectively.  For  accurate  representation  of  4// 
and  O',  large  number  of  mode  truncation  is  needed.  For  ex¬ 
ample,  for  the  RMS  error  of  the  toroidal  potential  (xtor)  to  be 
less  than  1%,  50  toroidal  modes  are  needed. 

Because  we  utilized  only  15  buoys  traced  the  Gulf  of  Mex¬ 
ico  circulation,  not  more  than  10-12  reconstruction  coeffi¬ 
cients  can  be  estimated  from  these  observations.  To  satisfy 
this  requirement  the  reconstruction  was  provided  on  186  Ju¬ 
lian  Day  of  model  calculations  because  at  this  time  we  can 
neglect  the  poloidal  current  and  approximate  the  toroidal  po¬ 
tential  through  10  toroidal  modes  with  the  accuracy  of  2.5%. 
The  second  reason  of  such  a  choice  is  that  the  reconstruction 
error  caused  by  heterogeneity  of  drifter  covering  of  the  area 
of  interest  was  considerably  reduced.  For  further  discussions 
see  Chu  et  al.  (2003a,  b). 

To  reconstruct  a  circulation  pattern  we  need  to  estimate 
mode  amplitudes  Ai,...,Aio  from  Lagrangian  data,  i.e. 
from  knowledge  of  drifter  velocity  (U,  V)  such  that 

||U-Vy(*  +  ¥')||2, 

+  1 1 V  +  Vt  (4'  +  4L)  1 1  >  min  (39) 

where  the  summation  in  the  Euclidian  norm  is  over  all  drifter 
positions.  The  reconstruction  coefficients  A i ,  . . . ,  Aio  were 
computed  by  the  Singular  Value  Decomposition  method 
(Engl,  et  al.,  1996). 

Figure  22  shows  the  model  circulation  at  186  Julian  Day, 
velocities  calculated  from  the  drifter  observations  for  the 
same  time,  the  reconstructed  currents  and  residual  circula¬ 
tion,  i.e.  difference  between  the  model  and  reconstructed  cur¬ 
rents.  It  is  noted  from  comparison  of  Figs.  22a  and  22c 
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Fig.  20.  RMS  errors  (xtor)  of  approximation  of  the  toroidal  circu¬ 
lation  through  1  (a),  3  (b).  6  (c),  10  (d),  30  (e)  and  50  (f)  toroidal 
EOFs. 


Fig.  21.  RMS  errors  (Xpol)  f°r  the  poloidal  circulation  approxi¬ 
mated  through  1  (a),  3  (b),  6  (c),  10  (d),  30  (e)  and  50  (f)  poloidal 
EOFs. 


that  the  spatial  structure  of  circulation  is  consistent  in  dif¬ 
ferent  scales  (the  Loop  Current  and  eddies)  between  the  re¬ 
constructed  and  model  fields.  The  residual  velocity,  i.e.  the 
reconstruction  error,  is  quite  high  in  the  central  part  of  the 
Gulf  of  Mexico  (Fig.  22d).  The  mean  relative  RMS  error 
between  the  model  and  reconstructed  fields  is  about  28%. 

7  Conclusions 

(1)  Lagrangian  predictability  of  the  high-resolution  Gulf  of 
Mexico  nowcast/forecast  system  developed  at  the  University 
of  Colorado  at  Boulder,  was  examined  using  55  Lagrangian 
drifters.  The  MPE  of  drifter  positions  and  velocities  was  es¬ 
timated  through  RMS  errors  and  irreversible-skill  time.  The 
approach  based  on  the  linear  tangent  model  was  found  un¬ 
satisfactory  because  in  mean  the  reference  period  can  not  be 
taken  longer  than  1  day.  Obviously  that  such  short  reference 
periods  are  of  no  interest  for  the  practical  applications. 

(2)  The  model  attractor  at  50  m  depth  was  analyzed 
through  EOF  decompositions  developed  for  toroidal  (non- 
divergent)  and  poloidal  (divergent)  components  of  circula¬ 
tion.  Lyapunov  dimensions  estimated  by  the  phenomenolog¬ 
ical  method  of  Syrovich  (1989)  held  the  following  inequali¬ 
ties 

15  <  £>toyrap  <16  1  <  Dpyap  <  2 

That  evidences  an  irregular  behavior  of  energetic  and  spec¬ 
tral  characteristics  of  the  circulation  and  probably,  exponen¬ 


tial  growth  of  prediction  error  caused  by  uncertainty  in  initial 
conditions  utilized  by  the  model. 

(3)  We  pointed  out  that  three  scenarios  are  possible  in  a 
model-drifter  comparison. 

(i)  The  long-term  (global)  predictability  is  referred  when 
the  model  reproduces  the  circulation  attractor  (the  ro¬ 
bust  regime)  in  small-scale  details  and  predicts  the  real 
drifter  trajectory  within  large  time  intervals. 

(ii)  In  the  short  and  intermediate  predictability  regime  the 
model  correctly  simulates  the  circulation  attractor  but 
not  all  small-scale  circulation  details.  The  drifter  trajec¬ 
tories  can  also  be  predicted  only  locally. 

(iii)  The  model  approximates  the  drifter  trajectories  only 
within  short  time  periods  (the  partial  predictability). 
Our  estimations  selected  the  second  scenario  as  one  re¬ 
alized  for  the  Gulf  of  Mexico  model. 

(4)  Statistics  of  MPE  is  nonGaussian  for  short,  intermedi¬ 
ate  and  long  times. The  MPE  grows  by  the  power  law  with 
exponent  equaled  to  3/2  in  the  range  12  h  <t<  20  days,  that 
corresponds  to  spatial  scales  between  20  km  and  300  km.  Ac¬ 
cordingly  to  terminology  of  the  theory  of  anomalous  diffu¬ 
sion  (Klafer  et  al.,  1987)  we  called  this  dynamical  regime  as 
“the  super-diffusion”  growth  of  prediction  error.  Obviously 
that  such  a  growth  indicates  the  fast  displacement  in  a  pair  of 
buoy-particle  in  areas  where  the  fixed  points  like  as  saddles 
may  be  identified  in  the  topology  structure  of  modeling  cir¬ 
culation.  In  our  opinion  the  power  law  of  model  prediction 
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Fig.  22.  Model  circulation  on  186  Julian  Day  (a),  velocities  calculated  from  the  drifters  for  the  same  time  (b).  the  reconstructed  circulation 
(c)  and  the  residual  circulation  (d). 


error  obtained  for  a  high-resolution  model  in  the  specific  ge¬ 
ographic  region  may  be  found  for  other  models  and  marine 
areas.  As  examples,  the  model  prediction  error  grows  along 
a  power  law  for  POP  model  of  the  North  Atlantic  (Ivanov  et 
al.,  2003). 

The  <7 — th  moment  of  MPE  is  scaled  with  power  exponent 
of  3 q  /4.  That  in  general  says  about  the  self-similar  feature  of 
model  prediction  skill.  The  PDF  of  MPE  holds  the  diffusion¬ 
like  Eq.  (34).  Universality  of  the  parameterization  of  Eq.  (34) 
to  describe  the  growth  of  prediction  error  for  other  high  res¬ 
olution  numerical  models  and  for  other  geographic  regions 
will  be  examined  in  a  separate  paper. 

(5)  The  PSP  clearly  shows  a  long  tail  stretched  into  large 
ITs.  This  tail  is  a  signature  of  long-term  correlations  between 
the  drifter  and  particle  trajectories  in  separate  individual  pairs 
of  drifter-particle.  Individual  Lagrangian  trajectories  can  be 
predicted  within  a  20-25  day  period  even  if  the  mean  La¬ 
grangian  prediction  time  does  not  exceed  4-5  days.  The 


best  fit  of  the  power  exponent  estimated  as  1.21  ±0.1 1  agrees 
very  well  with  the  theoretical  value  of  1 .33  obtained  from  the 
diffusion-like  MPE  model. 

(6)  The  obtained  results  has  the  following  practical  capa¬ 
bility.  First,  statistics  of  the  mean  prediciton  errors  can  be 
utilized  to  develop  numerical  scheme  for  Lagrangian  data  as¬ 
similation.  That  will  be  done  in  a  separate  paper. 

Second,  in  general  we  can  distinguish  between  different 
approaches  in  study  of  ocean:  Eulerian  (in  terms  of  veloc¬ 
ity  fields)  and  Lagrangian  (in  terms  of  trajectories  of  fluid 
particles).  Even  through  these  two  points  of  view  are  in  prin¬ 
ciple  equivalent,  the  relationship  between  predictabilities  as 
seen  in  the  realm  of  two  approches,  is  still  an  open  problem 
and  there  is  no  evidence  of  a  fixed  correspondence  between 
Eulerian  and  Lagrangian  chaotic  behaviours. 

It  was  pointed  out  in  Eremeev  et  al.  (1992)  that  if 
the  Lagrangian  correlation  scale  Ti  determined  from  the 
Lagrangian  data  exists  and  drifting  buoy  movements  are 
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described  as  pure  stochastic,  a  model  in  mean  may  not  pre¬ 
dict  Lagrangian  trajectories  for  time  intervals  (rpred)  longer 
than  this  scale,  i.e.  rpreci  <  7)  .  In  contrast,  if  drifter  move¬ 
ments  are  chaotic  (Lagrangian  chaos)  we  may  select  models 
allowing  to  predict  longer  that  Ti  (Kravtsov,  1989).  A  prac¬ 
tical  procedure  of  similar  selection  was  developed  by  Pires  et 
al.  (1996). 


Appendix  A  Computation  of  the  singular  values 


Let  us  introduce  the  matrix  Q—L  1  and  differentiate 
Eq.  (13)  with  respect  to  t: 


-  (  Q~l  ^  Q~l  sz,sz 


I  dSZ 

Q~x  - ,8 Z 

dt 


Q~l8Z,  )=0  (Al) 
dt 


d8Z 


Using  Eq.  (1 1)  Eq.  (A.l)  is  rewritten  as 


Q~l 


d  Q  , 


dQ 

or  ^  =  JiQ+QJ[ 

with  obvious  initial  conditions 
6  (to)  =  go 


Q-'SZ'SZ)  =0 


(A2) 

(A3) 

(A4) 


where,  go  is  a  non-zero  matrix  measuring  uncertainty  of 
model  initial  conditions. 

So,  for  the  given  reference  period  L— tp,  the  diagonal  ele¬ 
ments  of  the  matrix  g  determine  the  size  of  ellipsoid  of  pre¬ 
diction  uncertainty.  Non-diagonal  matrix  elements  identify  a 
position  of  this  ellipsoid  in  space.  Obviously,  the  geometrical 
center  of  the  ellipsoid  moves  along  the  reference  trajectory 
Zre f.  The  eigenvalues  and  eigenfunctions  of  the  matrix  g 
are  easily  recalculated  in  the  traditional  singular  values  and 
vectors  by  way  of  simple  algebraic  manipulations. 
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